function [W,param] = ecogTrainRegression(X,Y,regressionParam)

%[W,param] = ecogTrainRegression(X,Y,regressionParam)
%
% Purpose: Calculate the Multiple Input Multiple Output (MIMO) Wiener regression 
% matrix W that does Y=X*W, i.e. predicts Y (the movement or stimulus feature)
% from X (the ecog data). Ridge regression is used to estimate W.
%
% Y is regressed with multiple shifted version of X. Therefore, W extends in
% time and establishes a certain temporal relation between X and Y that can 
% be interpreted in terms of causality (see below). Heres X causal on Y means
% that effects in X have an influence on what happens later in time in Y. 
% 
% INPUT:
% X:        neuronal data, see ecogPrepareRegression
% Y:        The matrix of time series of external variables (measured movements,
%           or stimulus parameters). Variables change along the first
%           dimension and time increases along the second.
%
% regressionParam (optional): struct with fields
%           method - 0 for ridge regression (default) 
%                    1 for 'LARS-EN' 
%                    2 for LARS
%                    3 for LASSO (see
%                    http://www2.imm.dtu.dk/pubdb/views/publication_details.php?id=3897)
%                    normalization required!!
%                    4 for SVR (support vector regression, spider toolbox required) 
%           lambda (default = 0); 
%           epsilon (default = []); SVR epsilon
%           C (default = []); SVR C, set to inf for hard margin
%           normalize - 0 (default) no normalization
%                       1 zscore X 
%                       2 zscore X and Y
%           (Attention: if normalizing, keep param for prediction!)
%           
%
% OUTPUT:
% W:        The matrix W to predict Y from X in the equation Yhat=X*W.
%           The first enttry in W is an offset and must be discarde for
%           prediction.
% param: The zscore param for Y and X 
%           If regressionParam.normalize=true, use the param as input
%           for ecogWienerMIMORegressionPredict
%
% Requirements: 
% Currently the ststistics toolbox for ridge regression.
% LARS-EN,LARS,LASSO implememtation:
% http://www2.imm.dtu.dk/pubdb/views/publication_details.php?id=3897
% for elastic net computation

% 091129 JR wrote it based on Brian's code
% 2010: CR included LARS etc.
% 111111 CR seperated preparation of data and actaual train/prediction
% 120617 CR added SVR


%% Input check
if nargin <3,
    regressionParam=struct;
end
if ~isfield(regressionParam,'method'),
    regressionParam.method = 0;
end
if ~isfield(regressionParam,'lambda'),
    regressionParam.lambda = 0;
end
if ~isfield(regressionParam,'normalize'),
    regressionParam.normalize = 0;
end

% output check
if nargout >1,
    param.normalize=struct;
end

%% normalize, if enabled
if regressionParam.normalize>0,
    param.normalize.meanX = mean(X);
    param.normalize.stdX = std(X);
    X = (X-ones(size(X,1),1)*param.normalize.meanX)./(ones(size(X,1),1)*param.normalize.stdX);
end
if regressionParam.normalize>1,
    param.normalize.meanY = mean(Y);
    param.normalize.stdY = std(Y);
    Y = (Y-ones(size(Y,1),1)*param.normalize.meanY)./(ones(size(Y,1),1)*param.normalize.stdY);
end

%% Now the ridge regression
% We estimate the coefficients for each column in Y separately
if size(Y,1)~=size(X,1),
    error('Vector sizes of X and Y for regression doesnt match');
end


if regressionParam.method ==0,
    for k=1:size(Y,2)
        W(k,:)=ridge(Y(:,k),X,regressionParam.lambda,0);
    end
elseif regressionParam.method ==4, % support vector regression
        alg=svr;
        alg.optimizer='libsvm'; 
        if isfield(regressionParam,'epsilon') && ~isempty('epsilon'),
            alg.epsilon = regressionParam.epsilon;
        end
        if isfield(regressionParam,'C') && ~isempty('C'),
            alg.C = regressionParam.C;
        else
            alg.C = getDefC(X);
        end
        d=data(X,Y);
        alg.algorithm.verbosity=2;
        [res param.svr]=train(alg,d);
        W=[param.svr.b0 get_w(param.svr)];
else 
    if size(Y,2)>1,
        error('LARS-EN for multiple Y inputs currently not supported!');
    end
    if regressionParam.method ==1,
        %[s_opt, b_opt, res_mean, res_std] = crossvalidate(@larsen, 5, 100, X, Y(:,1),regressionParam.lambda);
        b1=larsen(X,Y(:,1),regressionParam.lambda);
    elseif regressionParam.method ==2,
        b1=lars(X,Y(:,1),'lars');
    elseif regressionParam.method ==3,
        b1=lars(X,Y(:,1),'lasso');        
    else
        error('Regression method unknown.');
    end
    %t1 = sum(abs(b1),2)/max(sum(abs(b1), 2));
    W=[zeros(size(b1,1),1),b1];
end
W=W';